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The phase diagram of prolate and oblate particles in the restricted orientations approximation 
(Zwanzig model) is calculated. Transitions to different inhomogeneous phases (smectic, columnar, 
oriented, or plastic solid) are studied through minimization of the fundamental measure functional 
(FMF) of hard parallelepipeds. The study of parallel hard cubes (PHC's) as a particular case is 
also included motivated by recent simulations of this system. As a result a rich phase behavior is 
obtained which include, apart from the usual liquid crystal phases, a very peculiar phase (called here 
discotic smectic) which was already found in the only existing simulation of the model, and which 
turns out to be stable because of the restrictions imposed on the orientations. The phase diagram 
is compared at a qualitative level with simulation results of other anisotropic particle systems. 
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I. INTRODUCTION 

Onsager first showed that the isotropic-nematic liquid 
crystal phase transition occurs in systems of anisotropic 
particles interacting via hard core repulsions . He stud- 
ied a system of hard spherocylinders in the limit of in- 
finite anisometry n = (L + D)/D — > oo (k is the sphe- 
rocylinder length to breath ratio) using the second virial 
form of the free energy, which in this limit is exact for the 
isotropic phase. The effect that higher virial coefficients 
have in the isotropic-nematic transition was later studied 
by Zwanzig, who introduced a model of hard prolate uni- 
axial parallelepipeds with axes oriented along the three 
perpendicular directions J2J. This peculiar model, which 
obviously treats the orientational degrees of freedom in 
an unrealistic way, has the advantage of being accessible 
to the calculation of higher virial coefficients up to sev- 
enth order in the infinite aspect ratio limit. He showed 
that including all these virials the isotropic-nematic tran- 
sition also occur, although the exact value of the coex- 
isting nematic density strongly depends on the order of 
the approximation. The Pade approximant generated by 
the truncated cluster expansion provides a much more 
stable sequence of the parameters which characterize the 
transition Q • This stability leaves little room for doubts 
regarding the existence of the transition in the model. 
The virial expansion resumed and expressed in the vari- 
able y = p/(l — pv), with p the number density of par- 
allelepipeds and v their volume, converges more rapidly 
than the traditional expansion in p, as was shown by Bar- 
boy and Gelbart for different hard particle geometries Q . 
Thus, the so-called 2/3 expansion of the Zwanzig model 
was applied to the study of the isotropic-nematic transi- 
tion as well as to the study of properties of its interface 
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[f| . For the latter the authors applied the smoothed den- 
sity approximation of the free energy functional in the 
spirit of Tarazona's weighted density approximation for 
the fluid of hard spheres [(| . 

The restricted orientations model for hard cylinders 
was also used to describe the structural properties of 
molecular fluids near hard walls or confined in a slit. 
This time the density functional was constructed from 
the bulk direct correlation function approximated by a 
linear combination of geometrical functions 0- 

Computer simulations of a variety of models of non- 
spherical hard core particles showed that the excluded 
volume effects could not only account for the stability of 
nematics, but also for the existence of liquid crystal inho- 
mogeneous phases such as the smectic 0] and columnar 
@ phases. Particularly the complete phase diagram of 
freely rotating hard spherocylinders |l0|, including not 
only smectic, but also a plastic solid phase and differ- 
ent oriented solid phases was calculated. Several den- 
sity functional theories, all of them based on weighted or 
modified weighted density approximations, are able to re- 
produce reasonably well the isotropic-smectic or nematic- 
smectic transitions ^2 113 m the whole range of as- 
pect ratios where the smectic is stable, and in some cases, 
transitions from the isotropic fluid to the plastic or ori- 
ented solid phases Q. In all these approximations the 
excess free energy is evaluated by integration of the free 
energy per particle of a reference fluid (typically spheres 
or hard parallel ellipsoids) evaluated at some weighted 
or effective density. In some cases, the employed weight 
is directly the normalized Mayer function between sphe- 
rocylinders 0, 0|; in others, it is calculated from the 
knowledge of the bulk correlation function of the refer- 
ence fluid For the latter case, the term proportional 
to the Mayer function enters into the integrand as a mul- 
tiplicative factor of the free energy per particle. The 
hard sphere free energy functional is recovered in both 
approaches as the limiting case of L = 0. 
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The fundamental measure theory (FMT) first devel- 
oped for hard spheres by Rosenfeld was another 
starting point for constructing a density functional for 
anisotropic particles. In its general formalism the ex- 
cess free energy density of the fluid is a function of some 
weighted densities obtained by convoluting the density 
profiles with weights which are characteristic functions of 
the geometry of a single particle whose integrals are the 
so-called fundamental measures: volume, surface area, 
and mean radius of the particles. Unfortunately, the 
Mayer function of two convex anisotropic bodies cannot 
be decomposed as a finite sum of convolutions of single 
particle weights ^(| , which is the keystone for construct- 
ing such a functional. Thus, the low density limit of the 
direct correlation function is no more the Mayer function. 

In spite of this, Chamoux and Perera have taken ad- 
vantage of Rosenfeld's extension of FMT to hard con- 
vex bodies by using it to compute the direct correla- 
tion function and patching out the low density limit with 
the exact Mayer function |17j . In this way they have 
obtained the equation of state for various convex hard 
bodies (such as hard ellipsoids, spherocylinders, and cut 
spheres), have predicted ordered phases and, recently, 
have study demixing in binary mixtures of rigidly ordered 
particles |l7llTl|. 

Following a similar procedure a density functional for 
anisotropic particles has been proposed which interpo- 
lates between the Rosenfeld's hard sphere functional and 
Onsager's functional for elongated rods. The resulting 
model was tested by calculating the isotropic-nematic 
transition in systems of hard spherocylinders and hard 
ellipsoids [T^|. 

Although the authors of this work suggest that the re- 
sulting theory can be applied to the study of inhomoge- 
neous systems, the huge computational efforts that their 
numerical implementations involve is the reason for the 
absence of any result in this direction. One way to cir- 
cumvent this difficulty is to reduce the continuous orien- 
tational degree of freedom to three discrete orientations 
(Zwanzig model). Implementing this idea some authors 
have recently applied the Zwanzig model to the study 
of interfacial properties of the hard rod fluid interacting 
with a hard wall or confined in a slit, for a one-component 
[20| and a polydisperse mixture [2l| , and also to the study 
of bulk and interfacial properties of hard platelet binary 
mixtures j2^. All these models are based on Onsager's 
density functional approximation. The increase of the 
number of allowed orientations in this functional partic- 
ularized for hard spherocylinders results in the presence 
of an artificial nematic-nematic transition in the one com- 
ponent fluid as the authors of Ref. |2^| have shown. This 
result indicates that certain cares must be taken in the 
direct extrapolation of the results obtained from this the- 
ory. 

FMF was also constructed for a mixture of parallel 
hard cubes combining Rosenfeld's original ideas with a 
dimensional cross over constraint |24j . The latter appears 
to be very important to describe correctly the structure 



of inhomogeneous fluids in situations of high confinement 
and to describe well the structural properties of the solid 
phase [2j| . The dimensional cross over has been used as 
an important ingredient to develop a density functional 
for a binary mixture of hard spheres and needles, assum- 
ing that the needles are too thin to interact with each 
other directly [2(j. 

Taking a ternary mixture of parallel hard cubes and 
scaling each species along one of the three Cartesian 
axes with the same scaling factor a FMF for the Zwanzig 
model is obtained. This functional has already been ap- 
plied to the study of the effect that polydispersity has 
on the stability of the biaxial phase in a binary mixture 
of rods and plates and on the relative stability of 
the smectic and columnar phases due to the presence of 
polydispersity [28| . 

The FMF for Zwanzig's model in the homogeneous 
limit coincide with the scaled particle theory and thus 
with the so-called 2/3 expansion which, as pointed out 
before, first began to be used in Ref. Q as a model 
to study the isotropic-nematic phase transitions in flu- 
ids of hard parallelepipeds. But this functional, through 
its minimization, also allows us to calculate inhomoge- 
neous density profiles. This functional has been applied 
recently to study the isotropic-nematic interface of a bi- 
nary mixture of hard platelets [29] . Its structural and 
thermodynamic properties resulting from the FMF min- 
imization show complete wetting by a second nematic. 
The same phenomenon was found in a binary mixture of 
hard spherocylinders [30| . 

The phase diagram for Zwanzig's model including the 
smectic, columnar, and solid phases has never been car- 
ried out, only spinodal instability boundaries have been 
traced [28]. The main purpose of this work is to obtain 
the complete phase diagram for this model and to com- 
pare the results with the only existing simulation of the 
lattice version of the model, which has been carried out 
for two different aspect ratios [3lJ . This will test the pre- 
dictive power of the FMF for anisotropic inhomogeneous 
phases. As a particular case, th e sy stem of parallel hard 
cubes will be studied. In Ref. [33 a bifurcation analy- 
sis and a Gaussian parametrization of the density profiles 
were used to calculate the free energy and pressure of the 
solid phase. Here a free minimization will be performed 
to calculate not only the solid but also the smectic and 
columnar phases and compare the obtained results with 
recent simulations of parallel hard cubes [13, 13 • 



II. FMF FOR ZWANZIG MODEL 

The FMF for hard parallelepipeds was already de- 
scribed in detail in Ref. j24[. A brief summary of the 
theory will be presented here putting emphasis on its 
numerical implementation to calculate the equilibrium 
inhomogeneous phases. 

A ternary mixture of hard parallelepipeds of cross sec- 
tion <7 2 and length L with their uniaxial axes pointing 
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to the x, y, or z directions is described in terms of their 
density profiles p^(r) (/j = x,y,z). Following the FMT 
for hard parallelepipeds in three dimensions the excess 
free energy density in reduced units can be written as 
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$ cxc (r) = $«(r) + $( 2 >(r) + $^(r), 
where the $( fc )'s are 
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i.e., they are sums of convolutions of the density profiles 
with the following weights: 



fe=i \ 

2fl(^-N) 



Ffc 



*(4-i* 



20 



Xj 



4°)(r), 



(6) 
(7) 

(8) 
(9) 



where the notation Xk (k — 1,2,3) for the x, y, and z 
coordinates has been employed. The functions 6(x) and 
8(x) are the usual delta Dirac and Hevisaide functions 
and aj t = a + (L — a)8^ with the Kronecker delta. 

The following constraints on the density profiles were 
imposed: (i) The solid phase has the simple parallelepi- 
pedic unit cell with uniaxial symmetry, i.e., the periods 
in the three spatial directions are d± for x, y and dn for 
z. The orientational director is selected parallel to z. (ii) 
The density profile of each species has the form 



/V(r) 



Plv 2^ a k II cos (ijkjXj) , (10) 

k=0 7 = 1 



where qj = 2vt / dj is the wave number along the j direc- 
tion, k = (fci, fc 2 , ks) is the vector defined by the recipro- 
cal lattice numbers, and n = (711,712,713) is the vector at 
which the harmonic expansion is truncated. Thus, Eq. 
HI U| l is the Fourier expansion of the density profiles ^(r) 
truncated at some n. This cutoff is selected in such a 
way that it guarantees small enough values of . The 



first Fourier amplitudes of all species are fixed to one 
(«o = !) and consequently J v ^ u drp^(v) = pj^ 
with V cc n = the unit cell volume, p the mean total 
density over the unit cell, and 7 M the occupancy proba- 
bility of species p in the unit cell, which obviously fulfills 
the condition X^7m = 

In the plastic solid phase these occupancy probabili- 
ties are 1 /3 for each species while they deviate from this 
value in the oriented solid phase. The uniaxial symmetry 
also implies that j x = j y = ~f±, j z = 711 = 1 — 27^ and 
p x (x,y,z) = p y (y,x,z), p z (x,y,z) = p z (y,x,z). Thus, 



the Fourier amplitudes verify a 



(x) 



and 



The total number of Fourier am- 



(*) _ («) 

plitudes [except the (0, 0, 0) term of all species] is reduced 
by these symmetries to N a = (n± + 1) (riy + l)(3nj_ + 
4)/2 — 2, (ni =rt2 = n.3 = nn) independent vari- 
ables. These variables together with j±, q± and qu span 
the variable space in which the FMF must be minimized. 

The density profiles of columnar and smectic phases 
are obtained from Eq. i|l(J|) substituting n = n±, 0) 
and n = (0,0, n\\). From the definition ©, Eqs. ©-© 
and the density profiles l|10(l . the weighted densities can 
be easily calculated resulting in 
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with v — La 2 the particle volume, (f>o(x) — cosx, 03 (x) = 
sinx/x, and £^ = qjkja^ /2. 

The substitution of Eqs. (|10|) and I jlljl into the free 
energy per unit cell 

[$id(r) + $cxc(r)], (16) 



t'cell Jv 



" ^ccll 

*id(r) = 5^p l ,(r)[ln(p |t (r)A3)-l], (17) 

with $id (r) the ideal part of the free energy density, and 
its minimization with respect to the N a + 3 variables 
allows the calculation of the equilibrium free energy and 
the density profiles of inhomogeneous phases. 
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To characterize the structure and orientational order 
of these phases the following total density and order pa- 
rameter profiles will be used: 



n! , _ , _ 3 [p a (r) +p y {r) 
Q{) 2 p(r) 



(18) 
(19) 



The selection of Q(r) as an order parameter is motivated 
by its uniaxial symmetry property Q(x, y, z) = Q(y, x, z) 
and its uniform limit value Q = 1 — 37j_ (—1/2 < Q < 
1), which coincides with the usual definition of the ne- 
matic order parameter for the Zwanzig model: Q = 
(7_L = 1/3) for the isotropic phase and Q — 1 ("fx = 0) 
for the perfectly aligned nematic phase. Although the 
solid and columnar phases might have local biaxiality 
[p x (x,y,z) ^ p y (x, y, z)], the integral over the unit cell 
of any previously defined biaxial order parameter is al- 
ways equal to zero as a consequence of the symmetries of 
the density profiles. 



III. PHASE DIAGRAMS 

The phase diagrams presented in this work were calcu- 
lated for a set of aspect ratios ranging from k = 0.1 
to k — 10, corresponding to the aspect ratios of the 
most anisotropic oblate and prolate parallelepipeds stud- 
ied here. The volume of all particles (cubes or prolate or 
oblate parallelepipeds) are fixed to 1 and thus the mean 
packing fraction rj is equal to the mean density p. From 
the equation v = La 2 = 1 the parallelepiped edge lengths 
L and a can be calculated as a function of the aspect ratio 
k = L/a as L = k 2 / 3 and a = k -1 / 3 . For each k, fixing 
the mean density p and using appropriate initial guesses 
for the N a + 3 variables with symmetries corresponding 
to the smectic, columnar or solid phases, the energy per 
unit cell l|16|) was minimized and thus the free energy 
for each phase was obtained. Varying p and repeating 
the former steps the free-energy branches of the differ- 
ent inhomogeneous phases have been calculated. The 
common tangent construction allowed the calculation of 
the coexisting densities between those phases in the case 
of first order transitions. To evaluate numerically the 
three dimensional integral of the free energy density i|ltj|) 
a Gauss-Chebyshev quadratures has been employed. 



A. Parallel hard cubes 

This subsection is devoted to the study of the parallel 
hard cube system (n — 1). The PHC equation of states 
of the fluid and solid phases as obtained from the FMT 
and the Monte Carlo simulation results are compared. 
While the solid phase is very well described with this 
formalism the exact location of the fluid-solid transition 
is very poorly estimated. The fundamental reasons of this 



difference are discussed here through a critical analysis 
of the fluid equation of state resulting from the FMF. 
It will be shown that possible modifications of the FMF 
slightly improve the location of the transition point at the 
expense of the correct description of the solid branch. 

In Ref. [27J the PHC fluid was already studied with 
the same FMF but using a Gaussian parametrization for 
the density profile. Through a minimization procedure 
and also from a bifurcation analysis a second-order fluid- 
solid transition was found at p = 0.3143 with a lattice 
period d = 1.3015 and a fraction of vacancies v = 0.3071 
|27j | . Recent simulations on the same system also showed 
a second-order transition to the solid but with very dif- 
ferent transition parameters p = 0.48 ± 0.02 in Ref. 33] 
and p = 0.533 ± 0.010 in Ref. 0. No evidence for the 
vacancies predicted by FMT was found, although the au- 
thors recognized that the vacancies might be suppressed 
by the boundary conditions in the small systems accessi- 
ble to simulations |34j . 

The main problem of the FMF for hard cubes is that 
it recovers in the homogeneous limit the scaled-particle 
equation of state, which overestimates the pressure 
calculated from the exact virial expansion up to seventh 
order. This expansion has a maximum at p w 0.6 and 
then goes down very quickly to reach negative values 
[35|. The poorly convergent character of the virial 
series makes it impossible to construct an equation of 
state for hard cubes, such as the Carnahan-Starling 
equation for the hard-sphere fluid, which estimates 
reasonably well all the known virial coefficients and 
diverges at close packing. On the other hand, it is 
well known that the FMF describes accurately the fluid 
structure in situations of high confinement, including 
the solid phase near close packing. For example, at 
high densities the functional recovers the cell theory, 
which is asymptotically exact when the packing fraction 
goes to 1, and also compares reasonably well with 
computer simulations |34j . These nice properties are a 
consequence of a fundamental restriction, namely, the 
dimensional cross-over |2 ll . imposed in the construction 
of the FMF. The latter implies that the functional in 
dimension D reduces to the functional in dimension 
D — 1 when the original density profile is constrained to 
D — 1 dimensions, i.e., p( D '(r) — p( D_1 ' (t)8(xd), where 
in is the coordinate that is eliminated on going from D 
to D — 1 dimensions. 

One possible procedure to improve the description of 
the uniform fluid of hard cubes at the level of the FMF 
is to follow the same method used in Refs. |36j and |37| . 
in which the hard-sphere Carnahan-Starling equation of 
state is imposed through the modification of the third 
term $( 3 ) [see Eq. of the excess free-energy density 
while keeping the exact density expansion of the direct 
correlation function up to first order. Unfortunately the 
absence of a good equation of state for the PHC fluid with 
the already mentioned properties makes this procedure 
less systematic compared to that of hard spheres [3(| |3?J • 
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Following this purpose the original excess free-energy 
density for hard cubes is now substituted by 

$c X c(r)=^/ fc (n 3 (r))$ (fc) (r), (20) 
fe=l 

with the //c(n3)'s selected in such a way as to keep the 
correct first order density expansion of the direct cor- 
relation function and to obtain the right virial expan- 
sion up to the seventh order of the PHC equation of 
state. As the original FMF for hard cubes gives the third 
virial coefficient correctly, these conditions imply that 
/i,2(n 3 ) - 1 + 0(^3) and fcim) ~ 1 + 0(n 3 ) for small 
713. Two further important conditions imposed on the 
/fe(«3)'s are their limiting behavior when the local pack- 
ing fraction tends to unity: lim„ 3 _>i / Q (n 3 ) = 1, which 
asymptotically guarantees the correct cell-theory limit, 
and the positive signs of their values, which guarantee the 
convexity of the fluid free energy. Unfortunately this pro- 
cedure breaks the dimensional cross-over property, but in 
principle should describe the fluid-solid transition in hard 
cubes better. 

Among all the functions f a 's that have been tried, even 
those which give better results [the particular case of 
/1, 2(713) = 1] are far from getting the transition point 
near the simulation one. In Fig. ^a) the scaled-particle 
equation of state, the improved equation of state 

^P = P + P^-$exc, (21) 

with <£> oxc being the uniform limit of Eq. (|20|l . and fi- 
nally, the symmetric Pade approximant of the seventh- 
order virial series are plotted. In the first two curves 
the bifurcation points are shown. The new bifurcation 
point calculated from Eq. i(27)|) is located at p = 0.3378, 
and the period and fraction of vacancies of the solid are 
d = 1.3249 and v = 0.2143. As can be seen from Fig. 
^a), the new equation of state still overestimates the 
fluid pressure, but to a lesser extent. Although the new 
functional gets a higher transition density and the frac- 
tion of vacancies decreases, there is still disagreement 
between theory and simulations. The equation of state 
of the PHC solid calculated from the minimization of 
the original FMF with respect to the Gaussian density 
profiles compare very well with simulations for densities 
p > 0.5, whereas the modified version underestimates the 
solid pressure. 

At this point the main conclusion that can be drawn 
is that the modification of the FMF in order to improve 
the description of the uniform fluid spoils the good de- 
scription of the solid phase. As the modification of the 
FMF was done at the expense of loosing the dimensional 
cross-over property (and this spoils the good description 
of highly inhomogencous phases), and the modified ver- 
sions do not show too many differences in the prediction 
of the fluid-solid transition, it is worthless to use them to 
study nonuniform phases. 
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FIG. 1: (a) Equations of state of the PHC liquid following the 
FMT (I/fmt), the modified version (Lmfmt), and the sym- 
metric Pade approximant. The circle and square represent 
the location of bifurcation points of the fluid-solid transition, 
(b) The equations of state of the solid phase from the original 
FMT (Sfmt) and from the modified version (Smfmt). The 
arrows represent the fluid-solid transitions predicted in Ref. 
pj^l (the lower value) and Ref. [34) (the higher value). Open 
and black diamonds are the simulation results from Ref. |33|| 
corresponding to the liquid and solid phases, respectively. 
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FIG. 2: Free energy per unit cell as a function of the mean 
density p. A linear function of p was subtracted from the free 
energy $*=<!> — mp — n to make clear the energy differences 
between the isotropic (J), columnar (C) and solid (S) phases. 



Setting q\\ = q± = q = 2n/d, 7 M = 1/3, and oc£' = /?k 
in the density profiles I|1U|I and minimizing the FMF, 
Eq. i|16|) . °f parallel hard cubes (n = 1) with respect to 
the Fourier amplitudes and the wave number q, the free 
energy per unit cell for solid [n — n(l,l,l)], columnar 
[n = n(l, 1,0)] and smectic [n = n(0, 0, 1)] phases were 
obtained. The results are shown in Fig. [2] From the 
isotropic liquid at the same density p — 0.3134 three in- 
homogeneous solutions: solid, columnar, and smectic, bi- 
furcate, with the solid phase being the stable one. While 



FIG. 3: Free energy per unit volume as a function of the mean 
density p. The involved phases are labeled as in Fig. [5] DSm: 
discotic smectic phase. The common-tangent constructions, 
which determine the coexisting densities labeled with different 
symbols, are also shown. 



the free energy difference between solid and columnar 
phases is relatively small, the smectic phase is clearly 
thermodynamically unfavorable. 

The number of Fourier amplitudes necessary to de- 
scribe adequately the density profile increases with the 
density, and thus the numerical calculations becomes 
more and more time consuming. Nevertheless, the sce- 
nario shown in Fig. [21 with the solid being the only 
stable phase, occurs at high densities as the simulations 
and cell-theory have confirmed [34| • The minimization of 
the FMF using a Gaussian parametrization of the density 
profiles of columnar and solid phases shows very similar 
quantitative results [34| ■ In fact the equation of state of 
the parallel hard-cube solid from FMT calculations with 
this parametrization compares very well at high densi- 
ties with simulations j2^. The results presented here 
are much more accurate than those obtained through the 
Gaussian parametrization. 



B. Prolate parallelepipeds 

This subsection is devoted to study the phase diagram 
of prolate particles (k > I). The results obtained from 
numerical minimization of the FMF of parallelepipeds 
with fixed k = 4.5 are shown in Fig. |3J The free en- 
ergies per unit volume of those phases which are stable 
in some range of densities are plotted. As can be seen 
the isotropic phase undergoes a first-order phase transi- 
tion to the so-called discotic smectic (DSm) phase. This 
peculiar phase is a layered phase (similar to the smectic 
phase) but with the long axes of the parallelepipeds ly- 
ing within the layers. There is no oricntational order in 
the layers, what means that the order parameter Q(z) 
reaches negative values at the positions of the density 
peaks. The density and order parameter profiles of the 
DSm phase at p — 0.3 are plotted in Fig. 21 The period 
in units of the small particle length is d/a = 1.2796 which 



FIG. 4: Density profile (solid line) and order parameter profile 
(dashed line) of the DSm phase at p — 0.3 (z* = z/d). 



means that the particles with long axes perpendicular to 
the layers (preferentially localized at the center of the 
interlayer space) intersect approximately three adjacent 
layers. 

Simulations of the Zwanzig model with k = 5 on a 
lattice showed an I-DSm transition at a density between 
0.47 and 0.55 .'if. Although the results were obtained 
for a lattice spacing of 1/3 (in units of the shortest par- 
ticle dimension) the simulations were repeated for val- 
ues 1/9 and 1/27 without changes in the stability of 
the DSm phase. Thus, the authors concluded that this 
layered phase may persist in the continuum limit |3l| . 
The difference in the transition density found from FMT 
(0.2868) and from simulations (~ 0.5) can be explained 
using two arguments: (i) As was already pointed out 
in Sec. [n] the FMF in the uniform density limit con- 
siderably overestimates the isotropic fluid pressure and 
thus the theory underestimate the transition densities 
between homogeneous and inhomogeneous phases, (ii) 
The transition densities should decrease upon decreasing 
the lattice spacing in simulations, as the results for the 
freezing of parallel hard cubes on a lattice (occurring at 
p = 0.568 for an edge length equal to two lattice spacings, 
at p = 0.402, for six lattice spacings, and at p = 0.314 
for the continuum) illustrate |38|. 

Increasing the mean density further the DSm phase 
undergoes a first-order transition to the columnar phase 
as Fig. shows. The restriction of parallelepiped ori- 
entations enhances the columnar phase stability even for 
prolate particles as a phase diagram, to be described be- 
low, will show. This phenomenon can be understood if 
the Zwanzig model is interpreted as a ternary mixture 
of particles. Simulations on a binary mixture of parallel 
spherocylinders with different aspect ratios (specifically 
2 and 2.9) show that, instead of a continuous nematic- 
smectic transition typical of the pure component sys- 
tem, the mixture exhibits a first-order nematic-columnar 
phase transition |39| . This result was explained by the 
poorer packing of rods of different lengths in the smectic 
phase as compared to that of rods of the same length. 
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FIG. 5: $* vs p for k = 4.5 (a) and k — 4 (b). OS: oriented 
solid phase. 



Simulations and theory show that one of the most im- 
portant effects that the aspect ratio polydispersity has 
on the phase behavior of hard spherocylinders l40l and 
binary mixtures of oblate and prolate particles |28l Elj 
is the enhancement of the columnar phase stability. All 
these results show that the columnar phase can be sta- 
ble even for mixture of particles with different shapes. 
Although the constituent particles of the Zwanzig model 
have the same shape, the restriction of their orientations 
changes strongly its relative packing and thus for some 
re's enhance the columnar phase stability with respect to 
other phases. 

At higher density the columnar phase exhibits a con- 
tinuous phase transition to an oriented solid phase of 
prolate parallelepipeds, as shown in Fig. EJa). The 
density and order parameter profiles of the columnar 
phase at the bifurcation point (p = 0.3748) are shown 
in Fig. [fj] The periods of the solid phase along the per- 
pendicular and parallel directions are dj_/a = 1.2690 
and du/L = 1.5170, respectively. From the equation 
p = (1 — (V ce n — d\dn being the unit cell vol- 

ume), the fraction of vacancies of the solid v can be calcu- 
lated as 0.0845. The continuous nature of the columnar- 
oriented solid transition changes to first order at some re 
between 4 and 4.5, as Fig. EJb) shows for re = 4. The 
order parameter Q(r) is very high in the unit cell ex- 
cept in its borders, where it exhibits small oscillations 
[see Fig. Efb)]. These oscillations are a consequence 
of the microsegregation of species "x" and "y" in the 
newly formed solid phase, which is preferentially formed 
by particles of species "z" localized around the position 
(x*,y*) = (0,0). This feature is shown in Fig. where 
the sum of the density profiles of species "x" and "y" 
[p_i_(r) = p x (r) + Py(r)] is plotted. While the columnar 
packing is responsible for the presence of the local max- 
ima at the center of the unit cell, the species "x" and 
"y" begin to segregate to the borders of the cell (±0.5, 0) 
and (0, ±0.5), respectively (see the four local maxima at 




FIG. 6: Density (a) and order parameter (b) profiles of the 
columnar phase at a density corresponding to the bifurcation 
point of Fig. EJa) (x* = x/d±; y* = y/d±). 



these positions) as the new solid phase is formed. The 
long axes of the perpendicular species lie on the lateral 
sides of the parallelepipedic unit cell, while their centers 
of mass are preferentially localized at the centers of these 
sides. 

The calculation of the free-energy branches for sev- 
eral stable inhomogeneous phases and the phase tran- 
sitions between them (as it was described for re = 4.5) 
has been carried out for 15 values of re (ten of them 
in the range 1 < re < 5 and five of them in the range 
5 < re < 10). The resulting phase diagram is plotted in 
Fig. |SJ The isotropic phase of prolate parallelepipeds 
with 1 < re < 3.5 undergoes a continuous phase transi- 
tion to the plastic solid phase. The transition points are 
joined with the spinodal line that has been calculated 
through the divergence of the structure factor. Notwith- 
standing that a functional minimization was carried out 
for each re to check the continuous nature of the transi- 
tions. The plastic solid is stable for re < 2.5 up to densi- 
ties around 0.5. At these values the numerical minimiza- 
tion turns out to be cumbersome because of the large 
number of Fourier amplitudes necessary to correctly de- 
scribe the inhomogeneous profiles. Thus, the high density 
part of the phase diagram (p > 0.5) has not been calcu- 
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FIG. 7: Sum of density profiles: p_i_(r) = pxijc) + p y (r) vs 
r* ± = (x*,y*) corresponding to the columnar phase shown in 

Fig. m 



lated with the numerical procedure described above. At 
higher densities a Gaussian-type parametrization of the 
density profiles is required, which obviously has a lower 
degree of accuracy. 

For re = 2.95 the plastic solid exhibits a very weak 
first-order phase transition to the discotic smectic phase 
(labeled as 1 in Fig. [SJ, and the latter a phase transi- 
tion to the columnar phase at higher densities. But the 
most representative region of the phase diagram where 
the discotic smectic is stable is for re around 4.5 where 
this layered phase exhibits a first-order phase transition 
to columnar phase (the shaded area of Fig. [8] limits the 
instability region against phase separation between both 
phases). For re between 4 and 5 the columnar phase un- 
dergoes a phase transition (first order for re = 4 and 
continuous for other values shown) to the oriented solid 
phase. The nematic phase begins to be stable for re > 5 
with its stability region bounded below by the first order 
isotropic-nematic transition and above by a continuous 
nematic-smectic transition. Finally, the smectic region is 
bounded above by a continuous transition to the oriented 
solid phase (see Fig. |HJ). 

Again the nematic-smectic transition points are joined 
with spinodal lines and for each re a minimization was 
carried out to check numerically the continuous charac- 
ter of the transition (the smectic solution begins to be 
stable right at the spinodal). In Ref. (53] a bifurca- 
tion analysis with the same functional was carried out 
to study the nature of the nematic-smectic transition. A 
thermodynamic and mechanical stability analysis showed 
that the nematic-smectic transition is first order, which 
is in contradiction with the numerical minimization re- 
sults presented here. A possible reason that justifies this 
contradiction could be that the iV-Sm transition is very 
weakly first order, so weak that the numerical accuracy 
used in the functional minimization can not decide about 
its nature. Another possibility is that the numerical ac- 
curacy failure is somewhere in the bifurcation analysis. 
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FIG. 8: Phase diagram of prolate parallelepipeds. Several 
phases are labeled like in Figs. |5J|^1 and[^](iV: nematic and 
PS: plastic solid) and the transition densities are labeled with 
different symbols (Star: isotropic; diamond: plastic solid; cir- 
cle: columnar; down triangle: discotic smectic; up triangle: 
smectic; and square: oriented solid). The shaded areas limit 
the regions of two phase coexistence. The transitions labeled 
by by 1,2, and 3 are first order in nature. 



A careful revision of this analysis is certainly called for 
in order to settle this point. 

The available simulation results for freely rotating hard 
spherocylinders show that the isotropic phase exhibits a 
transition to the solid phase for < re < 4.1 (the solid 
is plastic for re < 1.35 and oriented for 1.35 < re < 4.1) 
while the isotropic-smectic and nematic-smectic transi- 
tions begin at re = 4.1 and 4.7, respectively [Tjj [notice 
that for hard spherocylinders the length-to-breadth ratio 
is re = (L + D)/D]. We can see that, despite the differ- 
ent particle geometry and the restricted orientations of 
the Zwanzig model, the agreement for the threshold re at 
which spatial instabilities to the solid and smectic phase 
destabilize the homogeneous phases is rather good. Also 
the qualitative picture is similar: elongated rods form 
smectics, and more symmetric particles form solids. The 
main difference between them is that the Zwanzig phase 
diagram presents regions where the columnar and dis- 
cotic smectic phases are stable, a difference due to the 
restriction of orientations. 



C. Oblate parallelepipeds 

The phase diagram of oblate parallelepipeds (re < 1) is 
shown in Fig. The main differences after comparing 
the phase diagrams of prolate (Fig. |SJ) and oblate par- 
ticles are that in the latter: (i) The smectic is no more 
a stable phase, (ii) The region of columnar phase sta- 
bility is considerably larger, (iii) The stability region of 
the plastic solid is reduced (in fact this phase is stable 
only up to re -1 « 2.5) at the expense of that of the dis- 
cotic smectic phase, (iv) The transitions to the latter are 
strongly first order in nature (except for re -1 = 4.5). (v) 
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FIG. 9: Phase diagram of oblate parallelepipeds (the same 
symbols and labels of Fig. [S]are used). POS: Perfectly ori- 
ented solid. 



The oriented solid phase is replaced by a perfectly ori- 
ented solid in which "x" and "y" species are absent. This 
phase, after scaling in the z direction, is the same as the 
solid of parallel hard cubes. A solution from the FMF 
minimization with three dimensional spatial modulations 
and with 7_i_ 7^ has not been found in the parallelepi- 
pedic unit cell (the case of face-centered or body-centered 
cubic unit cells have not been tried here). 

Finally, similar by to what happens with prolate par- 
allelepipeds, the nematic phase begins to be stable at 
> 5. It undergoes a continuous phase transition to 
the columnar phase (the transition points of Fig. [5] are 
joined with the spinodal curve). 

The parallelepipeds with = 1.5 exhibit an inter- 
esting phase behavior. At low densities the isotropic 
phase destabilizes with respect to the columnar phase 
and not with respect to the PS phase. This example 
shows that the prediction for phase transitions using only 
the spinodal instability calculations can generate uncer- 
tainties about the possible symmetries of the inhomoge- 
neous phases. In fact these calculations do not allow to 
decide in this example if the new phase is a plastic solid 
or a columnar phase. Only by a complete minimization 
of the FMF, could it be concluded that the columnar 
phase is the stable one. 

In Fig. ^| (a) the density profiles of perpendicular 
[p^(z)] and parallel [p||(^)] species are shown for the dis- 
cotic smectic phase of oblate particles with k _1 = 2.5, 
while the order parameter profile is shown in Fig. HOt b). 
This discotic smectic phase coexists at p = 0.4244 with 
the plastic solid phase. The period in units of the large 
parallelepiped edge length is d/a = 1.2142. The random 
orientation of the uniaxial axes within the layers is con- 
firmed by the high negative values of the order parameter 
at the position of the density peak of the perpendicular 
species. The main difference between the DSm of Fig. ITU1 
and that of Fig. [7] is that the "z" species is now localized 
preferentially not at the center of the interlayer space, 




* 

z 



FIG. 10: Density (a) and order parameter (b) profiles of the 
DSm phase coexisting at p — 0.4244 with the PS. (a) Shows 
the density profiles of perpendicular (p±) and parallel (py) 
species. 



but near the smectic layers [see Fig. HOf a)]. exhibiting 
two local maxima at each side of the layer. This effect 
can be explained by the depletion force that the perpen- 
dicular species exerts on the parallel one. 

The iV-Sm (7Y - C) and the Sm-OS (C-POS) transi- 
tion lines of Figs. 00 and |§] converge asymptotically to 
p = 0.3143, the value of the fluid-solid bifurcation den- 
sity, as k — > 00 (k" 1 — > 00). The reason for this is that 
upon increasing k the number of rods (plates) with 

orientation perpendicular to the director becomes vanish- 
ing small, and then the system is, after rescaling the z 
direction, almost equivalent to a system of parallel cubes. 

Simulations of the Zwanzig model on a 15 x 15 x 15 
lattice with spacing 1/3 show that oblate parallelepipeds 
with dimensions 5x5x1 undergo a transition to a phase 
exhibiting a columnar structure |3ll | at a density some- 
where between 0.55 < p < 0.65. On increasing the sys- 
tem size to 30 x 30 x 30 the global columnar order dis- 
appears, but local correlations persist in the fluid with 
particle alignment distributed evenly among the three 
available orientations. In the same work the effect that 
orientational constraints have on the stability of the inho- 
mogeneous phases was studied. While a system of biaxial 
5x3x1 "tiles" without orientational constraints (except 
those inherent to the Zwanzig model) stabilizes in a smec- 
tic phase with the shortest edge lengths perpendicular to 
the layers, the system composed by "tiles" with all their 
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long edge lengths parallel to each other exhibits a phase 
transition to the smectic phase with these edge lengths 
perpendicular to the layers, similar to what is found here 
for uniaxial oblate parallelepipeds (the discotic smectic 
phase) . 

Simulations of hard cut spheres show that for n = 0.3 
there is an isotropic-solid transition, for k = 0.2 an 
isotropic-columnar transition (the isotropic phase might 
instead be a peculiar "cubatic" phase) and for n = 0.1 
a nematic-columnar one 0. From these results it can 
be concluded that the effect that the degree of particle 
anisotropy has on the symmetry of the stable phases for 
both cut spheres and hard parallelepipeds with restricted 
orientations, is qualitatively similar. 

IV. CONCLUSIONS 

The goal of this article has been the calculation of 
the phase diagram of the Zwanzig model for prolate 
and oblate parallelepipeds centering the attention on the 
phase transitions to inhomogeneous phases. For this pur- 
pose the fundamental measure functional for hard par- 
allelepipeds with restricted orientations has been used. 
This functional is exactly the same for any particle shape 
(prolate and oblate depending on k), which allows for a 
unified study of the phase behavior of both kinds of par- 
ticles. A free minimization of the functional was carried 
out with the only constraints of choosing a parallelepi- 
pedic unit cell and of imposing uniaxial symmetry in the 
inhomogeneous phases. The latter is justified by uniaxial 
symmetry of the particles. The degree of approximation 
to the exact density profiles was controlled by the cutoff 
imposed on the reciprocal lattice numbers in the Fourier 
expansion of the density profiles. 

A system of parallel hard cubes was separately studied, 
which was motivated by recent simulations on this sys- 
tem [33l 13 • Applying a modified versions of the FMF 
to improve the description of the PHC liquid, along the 
same lines as Refs. |36) and 37], the continuous transi- 
tion point to the solid phase and the equation of state 
of the solid were calculated from the divergence of the 
structure factor and from the functional minimization 
with respect to Gaussian density profiles. Although the 
transition density and fraction of vacancies change in the 
right direction, these results are still far from the simula- 
tions. In fact, the solid phase is poorly described by the 
new functional. The poor convergence of the PHC virial 
series does not make this procedure as effective as for 



hard spheres. Further refinement of the method and the 
proper inclusion of vacancies in simulations of the solid 
phase will probably improve the agreement between the- 
ory and simulations. 

The original FMF for PHC was minimized to study 
the relative stability of the smectic, columnar, and solid 
phases, starting at low densities from the bifurcation 
point. The solid phase is the only stable phase, followed 
by the columnar and the smectic (in order of energy sta- 
bility). At high densities the same behavior is shown 
from calculations using cell theory, functional minimiza- 
tion with Gaussian density profiles and computer simu- 
lations pi- 

The system of prolate and oblate parallelepipeds ex- 
hibits a very rich phase behavior. Apart from the plastic 
or oriented solid, smectic, and columnar phases, which 
are present also in systems of prolate (spherocylinders 
[Toj ') and oblate (cut spheres ||) particles, a new phase 
appears: the discotic smectic, the existence of which was 
confirmed by simulations |3l|]. The close relation be- 
tween the particle anisotropy and symmetry of the stable 
phases (elongated particles form smectics, flattened one 
form columnars and more isotropic particles form solids) 
which has been observed in simulations 0, 0, and 
experiments |43| is confirmed by this simple model. 

There are two important effects that the restriction 
of orientations has on the phase diagram topology: (i) 
The already pointed out stability of the discotic smectic. 
(ii) The stability of the columnar phase of prolate par- 
allelepipeds for some aspect ratios. The structural prop- 
erties of inhomogeneous phases that were found through 
functional minimization allow us to elucidate interesting 
effects such as the microsegregation behavior of differ- 
ent species in the solids and the depletion effect between 
particles in the smectics. Those findings endorse the pre- 
dictive power of the FMF in the description of highly 
inhomogeneous phases. 
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